/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Marcus Cuda
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
using System;

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    internal partial class DenseSvd : AbstractSvd
    {
        private double[] mMatrix;

        public DenseSvd(DenseMatrix matrix, bool computeVectors) : base(computeVectors)
        {
            mRows = matrix.Rows;
            mColumns = matrix.Columns;
            mMatrix = new double[mRows*mColumns];
            Buffer.BlockCopy(matrix.Data, 0, mMatrix, 0, mMatrix.Length*Constants.SizeOfDouble);
        }

        protected override void DoCompute()
        {
            int nm = System.Math.Min(mRows + 1, mColumns);
            mS = new DenseVector(nm);

            if (mComputeVectors)
            {
                mU = new DenseMatrix(mRows, mRows);
                mV = new DenseMatrix(mColumns, mColumns);
            }
            mConverged = Decompose();

            //adjust the size of s if row < columns.
            //we are using ported copy of linpack's svd code and it uses
            //a singular vector of length rows+1 when rows < columns.
            //the last element is not used and needs to be removed.
            //we should port lapack's svd routine to remove this problem.
            if (mRows < mColumns)
            {
                nm--;
                Vector tmp = new DenseVector(nm);
                for (int i = 0; i < nm; i++)
                {
                    tmp[i] = mS[i];
                }
                mS = tmp;
            }

            double eps = System.Math.Pow(2.0, -52.0);
            double tol = System.Math.Max(mRows, mColumns)*mS[0]*eps;
            mRank = 0;
            for (int h = 0; h < nm; h++)
            {
                if (mS[h] > tol)
                {
                    mRank++;
                }
            }

            //we no longer need the original matrix.
            mMatrix = null;
        }
    }
}